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The random phase approximation (RPA) for the electron correlation energy, combined with the 
exact-exchange (EX) energy, represents the state-of-the-art exchange-correlation functional within 
density- functional theory (DFT). However, the standard RPA practice - evaluating both the EX 
and the RPA correlation energies using Kohn-Sham (KS) orbitals from local or semi-local exchange- 
correlation functionals - leads to a systematic underbinding of molecules and solids. Here we 
demonstrate that this behavior can be corrected by adding a "single excitation" (SE) contribution, 
so far not included in the standard RPA scheme. A similar improvement can also be achieved by 
replacing the non-self-consistent EX total energy by the corresponding self-consistent Hartree-Fock 
total energy, while retaining the RPA correlation energy evaluated using KS orbitals. Both schemes 
achieve chemical accuracy for a standard benchmark set of non-covalent intermolecular interactions. 



In the quest for finding an "optimal" electronic struc- 
ture method, that combines accuracy and tractability 
with transferability across different chemical environ- 
ments and dimensionalities (e.g. molecules, wires/tubes, 
surfaces, solids), the treatment of exchange and correla- 
tion in terms of "exact-exchange plus correlation in the 
random-phase approximation (EX+cRPA)" 0, Q offers 
a promising avenue [3l-[l6j|. In this approach, part of the 
exact- exchange (EX) energy cancels exactly the spuri- 
ous self-interaction error present in the Hartree energy. 
The RPA correlation (cRPA) energy is fully non-local, 
whereby long-range van der Waals (vdW) interactions 
are included automatically and accurately [13] )• More- 
over, dynamical electronic screening is taken into account 
by summing up a sequence of "ring" diagrams to infi- 
nite order, making EX+cRPA applicable to small-gap or 
metallic systems where, for example, Hartrce-Fock (HF) 
plu s 2nd-order M0ller-Plesset (MP2) perturbation theory 
jp8j breaks down. 

The concept of cRPA dates back to the many-body 
treatment of the uniform electron gas in the 1950's 
0, Q ; and was later formulated [2(| within the context 
of density- functional theory (DFT) Recent years 

have witnessed a revived interest in EX+cRPA and its 
variants in quantum chemistry [ ja-lp , solid state physics 
10 12], and surface science I13H15I 1. Within the frame- 
work of Kohn-Sham (KS) DFT, EX+cRPA embodies 
an orbital-dependent functional that can in principle be 
solved self-consistently via the optimized effective poten- 
tial approach [2l|. This is however numerically very de- 
manding, and practical EX+cRPA calculations are com- 
monly performed in a post-processing fashion, where 
single-particle orbitals from a self-consistent DFT cal- 
culation in the local-density approximation (LDA), gen- 
eralized gradient approximations (GGAs), or alike, are 
used to evaluate both the EX and cRPA terms. Alterna- 
tively, one can formulate cRPA in terms of many-body 
perturbation theory (MBPT) based on a Hartree-Fock 



(HF) reference. 

Throughout this Letter we will adopt the following 
nomenclature: E F @SC is the total energy of the func- 
tional F, evaluated with the orbitals of a self-consistent 
(SC) scheme, e.g., HF, or the Perdew-Burke-Ernzerhof 
(PBE) [22[ GGA. The corresponding theoretical scheme 
is then labelled as F@SC. We also use the letter "x" or 
"c" in front of F or as a subscript of E F to refer to 
the exchange or correlation part of the scheme explic- 
itly. The functional F can be exact-exchange (EX), or 
additionally contain the RPA correlation (EX+cRPA), 
etc. For instance, £' EX @HF is the self-consistent Hartree- 
Fock energy, whereas the conventional RPA scheme based 
on PBE orbitals is referred to as (EX+cRPA)@PBE. 

The original (EX+cRPA) @PBE and (EX+cRPA) @HF 
schemes both exhibit systematic underbinding for a 
large variety of systems, including covalent molecules 
weakly bonded molecules J7 , [|| , solids [ll[ , and molecules 
adsorbed on surfaces 1 1 3l — tl 51 ) . Several attempts have been 



made to improve the accuracy of EX+cRPA. The earliest 
is the so-called RPA+ scheme [23| where a local correc- 
tion at the LDA/GGA level is added to cRPA. More re- 
cent attempts add second-order screened exchange (SO- 
SEX) 0,[24[) to make the entire approach self-correlation 
free, or invoke cRPA in a range-separated framework 
where only the long-range part of cRPA is incorporated 
Q. Among these, RPA+ improves total correlation 
energies considerably 25], but not binding energies 
The SOSEX correction performs well [t| [24| with con- 
siderable additional numerical effort. Range-separated 
RPA schemes also improve upon the standard EX+cRPA 
scheme 0, [H, [l6| , however, at the price of introducing em- 
pirical parameters in the approach. 

In this Letter, we offer a new perspective, based 
on MBPT, for going beyond cRPA, and show that a 
simple modification of the standard EX+cRPA scheme 
leads to a significant accuracy increase for molecular 
binding energies. We first illustrate our key idea us- 
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FIG. 1: (Color online) Panel (a): Binding energy curve 
for Av'2 computed with four RPA-based approaches, in 
comparison to the accurate reference curve by Tang 
and Toennies [2(J. Panel (b): Decomposition of the 
(EX+cRPA)@HF ((EX+cRPA)@PBE) binding energy of 
Ar 2 into individual contributions: EXOHF (EXQPBE) and 
cRPAQHF (cRPAOPBE). The difference between EX@HF 
and EX@PBE, and the SE@PBE term are also plotted. The 
vertical dashed line marks the equilibrium distance. Calcula- 
tions are done using FHI-aims [27l . [28j and Dunning's aug-cc- 
pV6Z basis [291 ]. The basis set superposition error (BSSE) is 
corrected here and in the following. 



ing the example of Ar 2 . The (EX+cRPA)@PBE and 
(EX+cRPA)@HF binding energy curves for Ar 2 are plot- 
ted in Fig. QJa). Both schemes show a significant under- 
binding behavior compared to the reference curve mod- 
eled by Tang and Toennies [26| based on experimental 
data. To gain more insight into the origin of the un- 
derbinding, the EX+cRPA binding energies are decom- 
posed into two contributions in Fig. [ljb): the exchange- 
only part and the remaining cRPA part. Inspection of 
the individual components reveals that ££ RPA @HF is 
(much) more repulsive than i^ RPA @PBE, whereas at 
the EX level £ EX @PBE is (much) more repulsive than 
£ EX @HF. The fact that ££ RPA @PBE is more attractive 
than E^ RPA @HF is easy to rationalize by inspecting the 
corresponding frequency-dependent polarizabilities. Ex- 
tensive benchmark calculations for 1225 molecular pairs 
30] show that asymptotic Cq dispersion coefficients de- 
rived from jE£ RPA @HF are systematically too small by 
approximately 40% [28J, while this error is only ~ 10% 
for E£ RPA @PBE. Adding AvdW corrections in an at- 
tempt to reduce the remaining error in cRPA@PBE (3lj 
only leads to minor changes in the binding energy at 
the equilibrium distance. What is more striking, how- 
ever, is the considerable difference in binding energies at 
the EX level — J B HF @HF- J B EX @PBE (plotted also in 
Fig. QJb) (red stars)). It amounts to ^6 meV at the 
equilibrium distance and is thus close to the deviation of 
the (EX+cRPA)@PBE binding energy from the reference 
value. 



From the viewpoint of Rayleigh-Schrodinger pertur- 
bation theory (RSPT), £ EX @HF and £ EX @PBE corre- 
spond to the sum of the zeroth and first-order terms in 
the perturbative expansions based on HF and PBE ref- 
erence state respectively 33|. The difference between 
£ EX @HF and £ EX @PBE must therefore be compen- 
sated by higher-order terms in the perturbation series 
since the final result should be independent of the refer- 
ence state, if all terms were summed up. The next term 

(2) 

in the series is the 2nd-order correlation energy Ec , to 
which only single and double excitation configurations 
contribute. Here we particularly examine the contribu- 

(2) 

tion of single excitations (SE) to £c , which can be ex- 
pressed [33j as 



E, 



a 



(1) 



to 



Here ipi and are the single particle orbitals and or- 
bital energies of the reference state, and / is the single- 
particle HF Hamiltonian - the Fock operator. A more 
detailed derivation of Eq. ([1]) is given in the supple- 
mentary material (where we simply follow the RSPT in- 
stead of the Gor ling-Levy PT [32|). As a consequence 
of the Brillouin theorem 33j, E^ E trivially vanishes 



for HF orbitals, but is in general non-zero for KS or- 
bitals. The contribution of E^ E evaluated with PBE 
orbitals (referred to as SE@PBE) to the binding en- 
ergy of Ar 2 is plotted in Fig. [TJb) (violet crosses). It 
amounts to 50% of the binding energy at the equilib- 
rium distance, and is close in magnitude to the contribu- 
tion from £ EX @HF-.E EX @PBE, and to the amount of 
under binding in the original (EX+cRPA)@PBE scheme. 
We therefore propose a new scheme by adding E to 
^ex+cRPA (subsequently referred to as EX+cRPA+SE). 
In Fig. [TJa) the resultant (EX+cRPA+SE)@PBE bind- 
ing energy curve is also plotted, which improves consid- 
erably over the (EX+cRPA)@PBE results, and is in close 
agreement with the Tang- Toennies reference curve. 

It appears that the quantitative agreement between 
£ C SE defined in Eq. Q and £ EX @HF-.E EX @PBE is a 
general feature. We found for a set of 50 atoms and 
molecules that the agreement typically ranges between 
70% and 100%, suggesting that replacing £ EX @PBE by 
_B EX @HF is an effective way to account for the SE con- 
tributions. This leads to a "hybrid-RPA" scheme, whose 
total energy is given by 



^hybrid-RPA _ £;EX ( 



!5HF + ^ RPA @PBE, 



(2) 



as an alternative to boost the accuracy of RPA. Fig.QJa) 
shows that the resultant binding energy curve is in almost 
perfect agreement with the reference curve. 

At this point, it is illustrative to take a closer look at 
the individual contributions to £ EX @HF-£ EX @PBE. In 
Fig. Owe further decompose the EX@HF and EX@PBE 
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FIG. 2: (Color online) Decomposition of the £ EX @HF and 
£ EX @PBE binding energies for Ar2 into their kinetic, elec- 
trostatic, and EX components. 



binding energies into their kinetic (T s ), electrostatic 
(_E clcc , external potential energy and Hartree energy com- 
bined), and EX components (E® x ) for Ar2. All three 
energy components behave quite differently for HF and 
PBE orbitals. The HF kinetic energy is purely repul- 
sive, whereas the PBE one exhibits spurious attraction 
at intermediate and large distances. The HF electro- 
static and exact-exchange energies, on the other hand, 
are purely attractive and decay to zero from below, while 
the corresponding PBE ones become repulsive in the in- 
termediate range and decay to zero from above at large 
distances. Since the PBE orbitals are less localized than 
their HF counterparts all three energy components de- 
cay much slower in PBE than in HF. The overall ef- 
fect is that i? EX @PBE becomes significantly more re- 
pulsive than i? EX @HF, resulting in the underbinding 
behavior of (EX+cRPA)@PBE. The more physical be- 
havior of EXOHF than EX@PBE at the EX level pro- 
vides a sound basis for the systematic improvement from 
(EX+cRPA)@PBE to hybrid-RPA. 

Indeed, the exceptional performance of the hybrid- 
RPA and (EX+cRPA+SE)@PBE schemes for rare-gas 
dimers carries over to many other molecular systems. As 
a second example we show results for the N2 molecule 
adsorbed on benzene (N2 ©benzene), which is an im- 
portant model system for studying molecular adsorp- 
tion on graphene and graphite surfaces |34| . We con- 
sider two possible configurations: N 2 placed parallel or 
perpendicular to the benzene plane. A successful the- 
oretical approach for this system must be able to de- 
scribe the delicate balance between electrostatic and 
dispersion interactions. We use FHI-aims [27], HU nu- 
meric atom-centered orbital basis (6s5p4d3/2<7 for C, 
O, N, and 5s3p2dlf for H ) augmented with gaussian 
diffuse functions from aug-cc-pV5Z to achieve conver- 
gence of the binding energy to within 1 meV. The re- 
sults shown in Fig. [3] are very similar to the rare-gas 
dimers: (EX+cRPA)@HF and (EX+cRPA)@PBE un- 
derbind significantly at the equilibrium distance, while 
hybrid-RPA and (EX+cRPA+SE)@PBE bring the bind- 
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FIG. 3: (Color online) Binding energies of the parallel and 
perpendicular configuration of N2 ©benzene as a function of 
the N2-C6H6 center-of-mass distance, calculated by different 
RPA-based approaches as well as MP2 compared to reference 
CCSD(T) calculations from Ref. 



TABLE I: Mean absolute error (in meV) and mean absolute 
percentage error (in parenthesis) of different RPA-based ap- 
proaches for the S22 databas e [3511 . CCSD(T) extrapolated to 
the complete basis set limit [3a] is taken as reference. 

H-bond Dispersion Mixed 

(EX+cRPA)@HF 45 ( 8.5%) 70 (43.9%) 34 (20.9%) 

(EX+cRPA)@PBE 57 (11.2%) 36 (21.8%) 24 (15.0%) 
(EX+cRPA+SE)@PBE 30 ( 6.0%) 18 (12.0%) 8 ( 5.5%) 
hybrid-RPA 18 ( 3.0%) 17 (10.0%) 8 ( 5.1%) 



ing energy into much closer agreement with the refer- 
ence curve computed with the coupled cluster method 
including single, double and perturbative triple excita- 
tions (CCSD(T)) [H]. In contrast, the traditional MP2 
method vastly overbinds the system. 

Finally we examine the performance of hybrid-RPA 
and (EX+cRPA+SE)@PBE for the S22 database of 
Jurecka et al. [351 ] . which represents a balanced bench- 
mark set for non-covalent interactions. The molecular 
dimers in this database can be divided into three groups 
of different bonding types: hydrogen-bonded, dispersion- 
bonded, and mixed complexes. We note that RPA in a 
range-separated framework has been applied to the S22 
database very recently . In Fig. [4] we plot the devia- 
tion from the CCSD(T) reference values [36( for the bind- 
ing energies in the S22 database [3f| for four RPA-based 
approaches and MP2. The basis set type and quality is 
the same as for N2@benzene. A detailed error analysis is 
presented in Table fl] 

We observe that the standard (EX+cRPA)@PBE 
scheme systematically undcrbinds all complexes. 
(EX+cRPA)@HF performs even worse for dispersion 
and mixed bonding, but better for hydrogen bonding. 
The latter case can be explained by the fact that the 
better performance of EX@HF dominates over the bad 
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molecules, one obvious deficiency of the 2nd-order SE as 
given by Eq. (p}, however, is that it is not well-behaved 
for systems with vanishing gaps. In such cases, we pro- 
pose to "renormalize" the SE contribution via a resum- 
mation of a geometrical series of higher-order diagrams 
involving single excitations (in the spirit of cRPA) . This 
leads to additional terms in the denorminator of Eq. (JTJ 
which prevent the possible divergence even when the KS 
gap closes. A brief derivation of this renormalized SE 
(RSE) scheme is presented in the supplementary mate- 
rial. Further details and benckmark calculations will be 
published elsewhere [28j |. 



FIG. 4: (Color online) Deviation from the CCSD(T) reference 
values [3(| for the binding energies of the S22 database [35j ] 
for RPA-based approaches as well as MP2. Positive errors 
correspond to overbinding and negative ones to underbinding. 



performance of cRPA@HF for hydrogen bonded systems. 
Again hybrid-RPA and (EX+cRPA+SE)@PBE correct 
the underbinding behavior of the standard EX+cRPA 
scheme, and improve the accuracy considerably. The 
hybrid-RPA scheme yields a mean absolute error (MAE) 
of 14 meV. The performance of (EX+cRPA+SE)@PBE 
is very similar to hybrid-RPA for dispersion and mixed 
bonding, albeit somewhat worse for hydrogen bonding. 
However, the mean absolute percentage error for hydro- 
gen bonding (6%) is still quite small. The accuracies 
achieved here compare favorably to the recently devel- 
oped vdW functional (vdW-DF) [37|, where the MAE for 
the PBE-based vdW-DF results for S22 [H is 54 meV. 
We also note that for covalent molecules the accuracies 
in the atomization energies are improved considerably 
by the two schemes. For instance, the MAE of the 
atomization energies of the G2-I set is reduced from 10.5 
kcal/mol to 6.2 kcal/mol by (EX+cRPA+SE)@PBE and 
6.3 kcal/mol by hybrid-RPA. 

To summarize, we have unraveled the origin of the un- 
derbinding that plagues the standard (EX+cRPA)@PBE 
scheme, which is mostly due to the too-repulsive 
nature of E EX @PBE rather than the (slightly) un- 
derestimation of the long-range dispersion force by 
££ RPA @PBE. This problem can be largely solved ei- 
ther by replacing S EX @PBE by the self-consistent HF 
energy i? EX @HF, or by adding a SE correction to 
the standard (EX+cRPA)@PBE approach. Particularly 
(EX+cRPA+SE)@PBE is a well-defined parameter-free 
scheme in which the SE term does not add any significant 
computational cost to the approach. In addition, the SE 
correction is compatible with other beyond- RPA schemes 
like RPA+ or SOSEX. We also like to emphasize that in 
both schemes the cRPA evaluated with KS orbitals is 
retained, which is essential for producing quantitatively 
correct asymptotics for vdW bonded systems. Despite 
its success for describing vdW and covalently bonded 
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DERIVATION OF THE SINGLE EXCITATION 
CONTRIBUTION TO THE 2ND ORDER 
CORRELATION ENERGY 

In this section we derive Eq. (1) that is presented 
in the main part of this Letter - the single excitation 
contribution to the 2nd-order correlation energy - from 
Rayleigh-Schrodinger perturbation theory (RSPT). The 
interacting TV-electron system at hand is governed by the 
Hamiltonian 



H = 



N r 
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where U e xt( r ) is a local, multiplicative external potential. 
In RSPT, H is partitioned into a non-interacting mean- 
field Hamiltonian H and an interacting perturbation H' , 



H = H° + H' 

N N 
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Here i) MF is some mean-field-type single-particle poten- 
tial, which can be non-local, as in the case of Hartree- 
Fock (HF) theory, or local, as in the case of Kohn-Sham 
(KS) theory. 

Suppose the solution of the single-particle Hamiltonian 
h is known 



(1) 



then the solution of the non-interacting many-body 
Hamiltonian H° follows 

H°\t> n )=EW\<I> n ). 

The |<&„) are Slater determinants formed from N of the 
spin orbitals \p) = \ip p ) determined in Eq. These 
Slater determinants can be distinguished according to 
their excitation level: the ground-state configuration 
$o), singly excited configurations |$"), doubly excited 
configurations |3>y), etc., where i,j, ... denotes occu- 
pied orbitals and a,b, . . . unoccupied ones. Following 
standard perturbation theory, the single-excitation (SE) 
contribution to the 2nd-order correlation energy is given 
by 
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where we have used the fact E^ — e\°^ = — e a . 

To proceed, the numerator of Eq. ^ needs to be 
evaluated. This can be most easily done using second- 
quantization formulation 
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where p, q, r, s are arbitrary spin-orbitals from Eq. ([T]), c\ 
and c q , etc. are the electron creation and annihilation op- 
erators, and (pq\rs) the two-electron Coulomb integrals 



(pq\rs) 



drdv 



r — r' 



The expectation value of the two-particle Coulomb oper- 
ator between the ground-state configuration $o and the 
single excitation $f is given by 

occ 

( $ o|r E^l rs )44 c ^l $ ?) = ^[{w\ap) ~ {ip\pa)] 



pqrs 



(3) 



where v HF is the HF single-particle potential. 

The expectation value of the mean-field single-particle 
operator v MF , on the other hand, is given by 



($0|E(^ MF |9>4 C *I $ "> = (Mv^a) 
pq 

Combining Eqs. ©, ®, and (gj), one gets 
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where Avt a is the matrix element of the difference be- 
tween the HF potential w HF and the single-particle mean- 
field potential i> MF in question. 

Observing that the -0's are eigenstates of h = 
— h V 2 + fcxt + v MF , and hence all non-diagonal elements 
{ipi\h°\ipa) are zero, one can alternatively express Eq. §5§ 
as 
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(6) 



where / is the single-particle HF Hamiltonian, or simply 
Fock operator. Thus Eq. (1) in the main paper is derived. 
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FIG. 1: Goldstone digrammatic representation of a selection 
of single excitation processes, (a): second order; (b-c): third- 
order; (d-g): fourth-order. Here vu = (i/'i|'C HF — v KS \ipi) and 

Vaa = <^|C HF -fi KS |^a). 
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For the HF reference state, i.e., when tr" = it 
the -0's are eigenstates of the Fock operator, and hence 
Eq. ([5]) (or ([6])) is zero. For any other reference state, 
e.g., the KS reference state, the i/Vs are no longer eigen- 
states of the Fock operator, and Eq. ([5]) is in general 
not zero. This gives rise to a finite SE contribution to 
the second-order correlation energy. In the following we 
assume v MF = v . 



a particle state a or a hole state i is, in the intermediate 
process, scattered into the same state by the perturba- 
tive potential Av = w HF — v KS . These terms dominate 
over the "non-diagonal" ones (2| and facilitate a simple 
mathematical treatment, because they form a geometri- 
cal series that can easily be summed up. Here we call the 
summation of this series renormalized single excitations 
(RSE). Following the textbook rule 0] for evaluating the 
Goldstone diagrams, one obtains 



E. 



RSE 



EE 

i a 
occ unoc> 

EE 

i a 
occ unoc> 

EE 



Aw*, 



3 \Av ta \ 2 (Av aa ~ Av u ) 



\Av ia \ 2 (Av aa - Av it ) 2 



'■.H.c uj..u ice | A ,2 



EE 



\Av ia \ 



e t - e a + Avu - Av a 



(7) 



The additonal term Ava — Av aa in the denominater of 
Eq. ([7|1, which appears to be negative definite in practical 
calculations, ensures that Ef- SE is finite and prevents pos- 
sible divergence problems even when the single-particle 
KS gap closes. In this context, we note in passing that 
a similar partial resummation of the so-called Epstein- 
Nesbet series Q of diagrams to "renormalize" the 2nd- 
order correlation energy has been performed by Jiang and 
Engel [Jj in the KS many-body perturbation theory. 



RENORMALIZED SINGLE EXCITATION (RSE) 
CORRECTION SCHEME 



ACKNOWLEDGEMENT 



The SE contribution as given by Eq. ([5]) corresponds to 
a correction to the correlation energy at 2nd-order, which 
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